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We report observation of 90° ferroelectric domain structures in transmission electron microscopy 
(TEM) of epitaxially-grown films of PbTiOs. Using molecular dynamics (MD) simulations based on 
first-principles effective Hamiltonian of bulk PbTiOs, we corroborate the occurance of such domains 
showing that it arises as metastable states only in cooling simulations (as the temperature is lowered) 
and establish characteristic stability of 90° domain structures in PbTiOs. In contrast, such domains 
do not manifest in similar simulations of BaTiOs. Through a detailed analysis based on energetics 
and comparison between PbTiOs and BaTiOs, we find that 90° domain structures are energetically 
favorable only in the former, and the origin of their stability lies in the polarization-strain coupling. 
Our analysis suggests that they may form in BaTiOs due to special boundary condition and/or 
defect-related inhomogeneities. 
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I. INTRODUCTION 



Ferroelectric transitions in perovskite oxides, such as 
BaTiOs, are fluctuation driven first-order phase transi- 
tions, and hence a state with spatially fluctuating order 
parameter can readily form as a result of certain me- 
chanical and electric boundary conditions^. A common 
example of such a state is the one with domains of fer- 
roelectric polarization with different symmetry equiva- 
lent orientations of order parameter that are separated 
by domain walls. Indeed, many properties of perovskite 
ferroelectrics depend on such domain structure, and it is 
being increasingly relevant at nano-scale 2 - - — . Naturally, 
the properties of a domain wall or an interface between 
adjacent ferroelectric domains depend on (a) symmetries 
and structural details of ferroelectric phases and (b) mi- 
croscopic couplings responsible for the ferroelectric phase 
transition. 

Perovskite oxides such as BaTiOs and PbTiOs are rep- 
resentative ferroelectric materials, although are quite dif- 
ferent from each other in terms of their phase transi- 
tions. While PbTiOs undergoes a single strongly first 
order phase transition from cubic to tetragonal struc- 
ture as temperature is lowered, BaTiOs exhibits a se- 
quence of three relatively weaker first-order phase transi- 
tions. A paraelectric phase of BaTiOs with cubic struc- 
ture transforms into a tetragonal ferroelectric phase at a 
Curie temperature, 393 K. Further cooling produces sud- 
den changes from a tetragonal phase to an orthorhom- 
bic phase at 278 K and from an orthorhombic phase 
to a rhombohedral phase at 203 K^i. On the other 



hand, PbTiOs exhibits a phase transition from a para- 
electric cubic phase to a ferroelectric tetragonal phase at 
Tq = 763 K and remains tetragonal down to K— ~— . 

While the domain structures in ferroelectrics have been 
revealed experimentally 7 -*^ - — , detailed in situ experi- 
mental analysis of the domains and domain walls is quite 
challenging and the temperature dependence of dynamics 
of domain structure is not well understood. In the case 
of PbTiOs, experimental studies of domains are further 
more difficult as the sample needs to be heated over its 
high transition temperature Tq — 763 K and such heating 
leads to evaporation of Pb ions changing the composition 
of the sample^. 

Ferroelectric phase transitions in perovskite oxides 
in bulk and thin film have been investigated by com- 
puter simulations such as phase- field method 27 , Monte 
Carlo simulations 2 ^, and molecular dynamics (MD) 
simulations—. Recently, Nishimatsu et al. have devel- 
oped a fast and versatile MD simulator of ferroelectrics 
based on first-principles effective Hamiltonian22 which 
can be used in systematic studies of bulk as well as thin 
films. They have studied BaTiOs bulk and thin-film ca- 
pacitors and obtained results showing good agreement 
with experiments. Their MD simulations of BaTiOs un- 
der periodic boundary condition (PBC) for bulk did not 
show any domain structures, as there is no depolariza- 
tion field in the PBC of bulk. Simulations of thin-film of 
BaTiOs only show 180° domain structures^,, though 90° 
domain structures are widely seen in experiment o 1 1 ' 13 . 
One of the advantages of MD simulations compared 
to Monte Carlo simulations is its ability to simulate 



time-dependent dynamical phenomena, e.g. MD simu- 
lation can be used to study the evolution of ferroelec- 
tric domains as a function of time during heating-up and 
cooling-down simulations. 

In this paper, we report heating-up and cooling-down 
molecular-dynamics (MD) simulations of bulk PbTiC>3 
to understand our observation of 90° domain structures 
in epitaxially-grown sample of PbTi03. In Sec. UH we 
present experimental details for the sample preparation 
and we show a transmission electron microscope (TEM) 
image of 90° domain structure in PbTi03 film. We briefly 
explain the first-principles effective Hamiltonian and de- 
tails of MD simulations in Sec. IIIII and we present our 
results and analysis of heating-up and cooling-down MD 
simulations in Sec. IIVI We finally summarize our work 
and conclusions in Sec. [Vj 



II. EXPERIMENTAL DETAILS AND 
OBSERVATIONS 

A. Sample preparation and Methods of TEM 

An epitaxial PbTiC>3 thick film, with film thicknesses 
of about 1200 nm, was grown on the SrRuOs/SrTiOs 
substrate at 873 K by pulsed metal organic chemical va- 
por deposition (pulsed-MOCVD) method. SrRu03 was 
deposited on (100) SrTiOs by rf-magnetron sputtering 
method. The detail of film preparation technique is de- 
scribed elsewhere2ii22. The TEM specimens were pre- 
pared with focused ion beam (FIB) micro-sampling tech- 
nique. Damage layers, introduced during FIB microfab- 
rication, were removed by low-energy Ar ion milling at 
0.3 kV. JEM-2000EXII was used for TEM observations. 
TEM observations were performed at room temperature. 

B. Observed TEM image 

In Fig. [TJ we show a bright field TEM image of a 
PbTiC>3 thick film, taken with the incident electron beam 
parallel to the [100] axis of the PbTiC>3. The inset in the 
image is the corresponding selected-area electron diffrac- 
tion pattern. This bright-field TEM image is taken with 
the scattering vector g n excited. Here, a- and c-domains 
are those with polarization along a and c axes of the 
PbTiC>3 parallel and perpendicular to the substrate, re- 
spectively. This domain configuration is typical and very 
commonly seen for tetragonal PbTiOs films. The do- 
main size is about 50-200 nm. Such 90° domain configu- 
rations, similar to that in Fig. [I] have been also observed 
in BaTiO ^ 11 ' 13 . Boundaries between a- and c-domains 
are seen as black lines, as the TEM sample is slightly 
tilted. From this image, we can not discuss the width of 
domain walls between a- and c-domain. High-resolution 
TEM observation has revealed that the width of 90° do- 
main walls is 1.0 ±0.3 mn2. To this end, high-resolution 
TEM observation was conducted in order to reveal the 
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FIG. 1: A bright-field TEM image of a PbTiOs thick film 
taken with an electron incident parallel to the [100] of 
PbTiOs. Scripts a and c in this figure denote a-domain and 
c-domain, respectively. 



width of domain walls between a- and c-domain. Fig. [5] 
shows the high- resolution TEM image of PbTiC>3 film, in- 
dicating that the width of domain walls is corresponding 
to 1 or 2 unit cells. 

Before our computational study, it should be worth 
mentioning that 90° domains have been often observed 
in both BaTiC>3 and PbTiOs, and are both ferroelectric 
and ferroelastic in nature. Typically, 90° domains are 
formed in epitaxial ferroelectric and ferroelastic films in 
order to relax the strain resulting from lattice mismatch 
with the substrate at and below Tc i 34 ' 35 They nucleate 
at misfit dislocations formed above Tq. Their growth is 
accompanied with the introduction of the additional dis- 
location perpendicular to the misfit dislocations and the 
dissociation of the dislocations into two pairs of partial 
dislocations around an anti-phase boundary 23 . 



III. MOLECULAR DYNAMICS SIMULATIONS 

A. Effective Hamiltonian 

Heating-up and cooling-down molecular-dynamics 
(MD) simulations are performed using first-principles ef- 
fective Hamiltonian^^i 3 ^^, 
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FIG. 2: A high-resolution TEM image of a PbTi0 3 thick 
film taken with an electron incident parallel to the [100] of 
PbTiOa. The 90° domain boundary is shown by the chain 
line. 



M * M* 

HCS = ^ole £ . 1{R) + ^^Y^WKR) 

R,a R,a 

+ V scli ({u}) + V dp \{u}) + V shmt ({u}) 
+ y 01 ^ 1101110 ^, . . .,r? 6 ) + y clas ' inho ({ti;}) 
+ V cmp > homo ({u}, m ,..., Tfe) + V co ^ inho ({u}, {w}) , 

(1) 

where u = u(R) and w = w(R) are, respectively, the 
local dipolar displacement vector and the local acoustic 
displacement vector of the unit cell at R in a simulation 
supercell. a(= x, y 7 z) is the Cartesian directions. Braces 
{} denote a set of u or w in the supercell. r]i,...,rie 
are the homogeneous strain components, -^dipoic an< ^ 
M* coustic are the effective masses for u and w, therefore, 
first two terms in Eq. (JIJ are kinetic energies of them. 

^/sclf ^/dpl ^short ^clas, homo ^/clas. inho ^coup, homo 

and l/ cou P> lnho are a local-mode self-energy, a long-range 
dipole-dipole interaction, a short-range interaction, a ho- 
mogeneous elastic energy, an inhomogeneous elastic en- 
ergy, a couping between {«} and r)i, . . ., rje, and a couping 
between {u} and {w}, respectively. More detailed ex- 
planation of symbols in the effective Hamiltonian can be 
found in Refs. Ho| and H3- We take all the parameters of 
the first-principles effective Hamiltonian for PbTiOa from 
the earlier work 3 ^. However, the form of the effective 
Hamiltonian we use in our simulation o 30 ' 37 is slightly dif- 
ferent from the one used to get the parameters in Ref. l38l . 



TABLE I: The parameters of first-principles effective Hamil- 
tonian for PbTiC>3 used in our simulations are given in the 
second column and how these parameters are related with 
the parameters in the previous work2& are shown in the third 
column. 



parameters 


value 


relation 


a [A] 


3.969 


a 


Bii [eV] 


117.9 


C11 


B 12 [eV] 


51.6 


C12 


B 44 [eV] 


137.0 


C44 


B lxx [eV/A 2 ] 


-114.02 


2(5o+3o) 


B lyv [eV/A 2 ] 


-13.67 


250 


B iyz [eV/A 2 ] 


-22.67 


92 


a [eV/A 4 ] 


27.83 


B + C 


7 [eV/A 4 ] 


-34.48 


-2C 


fei [eV/A 6 ] 


-42.42 


D 


k 2 [eV/A 6 ] 







fca [eV/A 6 ] 







k 4 [eV/A 8 ] 


156.43 


E 


m* [amu] 


100.0 




Z* [e] 


10.02 


Z* 




8.24 




K 2 [eV/A 2 ] 


1.170 


A 


31 [eV/A 2 ] 


-1.355 


2a T 


32 [eV/A 2 ] 


4.986 


2a L 


h [eV/A 2 ] 


0.222 


bi + b a 


34 [eV/A 2 ] 


-0.018 


2b t2 


h [eV/A 2 ] 


0.398 


bi - b t i 


36 [eV/A 2 ] 


-0.083 


2(c,+2c t ) 
3 


37 [eV/A 2 ] 


-0.204 


2( c ,-2c t ) 
?■ 



The new parameters can be easily derived from the pre- 
vious ones. We list values of all the parameters used in 
our simulations and how they are related to the previous 
work££ in Table HI 



B. Simulation Details 

Heating-up and cooling-down MD simulations 
are performed with our original MD code feram 
(http://loto.sourceforge.net/feram/). Details of 
the code can be found in Ref. 30. Temperature is 
kept constant in each temperature step in the canon- 
ical ensemble using the Nose-Poincare thermostat. 39 
This simplectic thermostat is so efficient that we 
can set the time step to At — 2 fs. In our present 
MD simulations, we thermalize the system for 20,000 
time steps, after which we average the properties for 
20,000 time steps. We used a supercell of system size 
N = L x x L y x L z = 32 x 32 x 32 and small temperature 
steps in heating-up (+1 K/step) and cooling-down 
( — 1 K/step) simulations. The heating-up simulation 
from 100 K to 900 K is started from an z-polarized initial 
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FIG. 3: (Color online) Simulated temperature dependence of 
lattice constants of PbTi03 in heating- up (red solid lines) and 
cooling-down (blue dashed lines) molecular-dynamics simula- 
tions. 



configuration generated randomly: (u x ) — (u y ) = 0, 
(u z ) = 0.33 A, (ul) - {u x f = (ul) - {u y f = (0.045 A) 2 , 
and (u 2 z ) - (u z ) 2 = (0.021 A) 2 , where brackets denote R- 
average in supercell (u a ) — jjJ2 R u a (R)- The cooling- 
down one from 900 K to 100 K is started from random 
paraelectric initial configuration: (u x ) = (u y ) = (u z ) = 
and (u 2 a ) - (u a ) 2 = (0.15 A) 2 . 



IV. RESULTS AND DISCUSSION 

From the temperature dependence of averaged lat- 
tice constants (shown in Fig. [3]), a tetragonal-to-cubic 
ferroelectric-to-paraelectric phase transition is clearly ob- 
served in the heating-up simulation at 677 K. However, 
a strange behavior in lattice constants is found in the 
cooling-down simulation at T=592 K. The average tem- 
perature of these two transition temperatures (634 K) 
is in good agreement with the earlier Monte Carlo sim- 
ulations and slightly lower than the experimental value 
Tq = 763 K. Indeed, the observation of an orthorhombic 
phase during cooling-down simulations is intriguing. 

To understand this interesting behaviour of lattice con- 
stants in cooling-down simulations, we perform a detailed 
analysis of the configurations (snapshots) during our MD 
simulations. From a snapshot of dipoles in the supercell 
(shown in Fig.^J), we find that the apparently orthorhom- 
bic nature of the phase is due to a 90° domain struc- 
ture. Although the 4 unit cell = 1.6 nm of the domain 
size is much smaller than experimentally observed ones 
as shown in Sec. Ill B[ the width of a simulated domain 
wall estimated to be ~ 1 unit cell is in good agreement 
with our experiment. Each domain has the a = b < c 
tetragonal structure of PbTiC>3, but their average value 
in whole crystal gives smaller d than c and larger a! than 
a. The lattice constant b' has almost the same values as 



a, because the polar directions of two kind of domains 
are perpendicular to the 6'-axis. It should be noted that 
this domain structure is found in MD simulations in bulk 
under periodic boundary condition (PBC), but we have 
not simulated thin films. Under the PBC, there is no 
depolarization field inside the bulk. Moreover, this do- 
main structure can be easily reproduced in cooling-down 
simulations from any random paraelectric initial config- 
urations and any seeds for the pseudo random number 
generator—. There was no evidence for such unusual be- 
havior in simulations of bulk BaTiO ^ 30 ' 37 . 

To understand the reason of stability of the 90° domain 
structure seen here, even in bulk PbTiC>3, we compare 
"total energy surfaces" between single and 90° domain 
structures. The total energy surface of single domain 
structure with [001] polarization is the same as in Refs.l36l 
and HO- For the total energy surface of 90° domain struc- 
ture, we focus on a snapshot of the supercell at 300K 
shown in Fig.@]and represent it with {m9o°30ok(-R)}- We 
now obtain a sequence of configurations by multiplying 
a factor ,. A " for all R 



(")o 



u(R) 



90°,300K 



" M 90° 300K(-R) , 



(2) 



and compute total energy as a function of u, where 
(w)go° 300K is the averaged length of dipoles in the 300 K 
snapshot 

(m)90°,300K = ^2 l M 90 o ,300K(J?)| = 0.32 A . (3) 
ft 

Calculated total energy surfaces of single and 90° do- 
main structures for PbTiC>3 are shown in Fig. [3] with solid 
and dashed lines, respectively. For comparison, those for 
BaTiC>3 are also plotted assuming the same 90° domain 
structure by using the set of parameters in H eS listed in 
Ref. [13. While the 90° domain structure of PbTi0 3 ex- 
hibits a minimum at u ^ 0, that of BaTi03 costs energy. 
This is why the 90° domain structures can be found in 
simple cooling-down simulations of PbTiC>3, but not in 
those of BaTi0 3 . 

Minima are indicated with (a)-(c) in Fig. [5j To un- 
cover the origin of this contrasting behaviour, interac- 
tion energy terms at the minimums are listed in TA- 
BLE [Til For the BaTi0 3 , because there is no minimum, 
interaction energies of configuration at (d) in Fig. [5] (of 
u = 0.10 A) are listed. 

From TABLE [UJ It is clear that the energy losses in 

yharmonic^ ^clas.inho^ and ycoup.homo j n f ormm g tne 90° 

domain structure in PbTi03 are compensated by the en- 
ergy gains in ^/ un harmonic ^clas, homo \T cou P> hiho 

contrast, such recovery is not sufficient in BaTi03 to form 
an / minimum or non trivial 90° domain structures. 
In stabilization of the 90° domain structures in PbTiC>3, 
our analysis conclusively highlights the role of two mi- 
croscopic interactions: (a) lower elastic energy cost aris- 
ing from the smaller strain from compensation along c 
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FIG. 4: (Color online) Snapshot of three "sides" of the 32 x 32 x 32 supercell at T — 300K in a cooling-down simulation of 
PbTi03. Dipole moments of each site are projected onto the side planes and indicated with arrows. Dipoles of u z > 0.2A are 
indicated with red color. Dipoles of u z < 0.2 A are indicated with green. Crystalline directions are indicated with a', &', and c' 
as indicated in Fig. [3] 



TABLE II: Comparison of interaction energies V % in the effective Hamiltonian of Eq. fT]) for two kinds of domain states of 
PbTiOa and BaTi0 3 . y harmonic is the sum of V dpl , V Bholt , and the harmonic terms in V sclf . v unh * lmonic is the unharmonic 
terms in V . Unit of energy is eV. 



x of 

interaction 
energy V x 


PbTiOs PbTiOa 
(a) single (b) 90° 
domain domain AV* 


BaTi0 3 BaTi0 3 
(c) single (d) 90° 
domain domain AV l 


harmonic 

unharmonic 

elas,homo 

elas,inho 

coup, homo 

coup,inho 

total 


-0.22106 -0.14393 +0.07713 
0.33435 0.17249 -0.16186 
0.21946 0.04441 -0.17505 
0.00000 0.05171 +0.05171 

-0.43891 -0.08881 +0.35010 
0.00000 -0.10342 -0.10342 

-0.10616 -0.06755 +0.03861 


-0.03886 0.03541 +0.07427 
0.04754 0.00695 -0.04059 
0.02609 0.00131 -0.02478 
0.00000 -0.00262 -0.00262 

-0.05218 0.00089 +0.05307 
0.00000 -0.00178 -0.00178 

-0.01741 0.04016 +0.05757 


u[A] 


0.34 0.29 


0.16 0.10 
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FIG. 5: (Color online) Calculated total energy surfaces 
for PbTi0 3 single domain (solid line), PbTi0 3 90° do- 
main (dashed line), BaTiOs single domain (dotted line), and 
BaTiOs 90° domain (chain line). 



V. SUMMARY 



In this article, TEM observation of PbTi03 thick films 
revealed 90° domain structures, which have been often 
observed in BaTiOs. The domain size perpendicular to 
domain boundaries was 50-200 ran. The width of domain 
wall was corresponding to 1 or 2 unit cell. 

We also have performed heating-up and cooling-down 
MD simulations of PbTiOs. In cooling-down simulation, 
90° domain structure is found to form spontaneously. 
By comparing "total energy surfaces" of single and 90° 
domain structures, we understand that a 90° domain 
structure is metastable in bulk PbTiOs, but not in bulk 
BaTiOs. The origin of this contrast is traced to sig- 
nificantly larger polarization-stran coupling in PbTiOs. 
Hence, while 90° domain structures can form sponta- 
neously in PbTiOs, they seem to arise in BaTiOs mostly 
from special boundary conditions and/or defect-related 
inhomogeneities. 
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